Draft version April 17, 2009 

Preprint typeset using I^TgX style emulateapj v. 9/08/03 



A NEW SECULAR INSTABILITY OF ECCENTRIC STELLAR DISKS AROUND SUPERMASSIVE BLACK HOLES, 

WITH APPLICATION TO THE GALACTIC CENTER 

Ann-Marie Madigan, Yuri Levin and Clovis Hopman 

Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden 
Draft version April 17, 2009 



o 
o 

< 



O 



> 

in 
o\ 
m 
m 

oo 
O 



X 



Abstract 

We identify a new secular instability of eccentric stellar disks around supermassive black holes. We show that 
retrograde precession of the stellar orbits, due to the presence of a stellar cusp, induces coherent torques that 
amplify deviations of individual orbital eccentricities from the average, and thus drive all eccentricities away 
from their initial value. We investigate the instability using A^-body simulations, and show that it can propel 
individual orbital eccentricities to significantly higher or lower values on the order of a precession time-scale. 
This physics is relevant for the Galactic center, where massive stars are likely to form in eccentric disks around 
the SgrA* black hole. We show that the observed bimodal eccentricity distribution of disk stars in the Galactic 
center is in good agreement with the distribution resulting from the eccentricity instability and demonstrate 
how the dynamical evolution of such a disk results in several of its stars acquiring high (1 — e ^ 0.1) orbital 
eccentricity. Binary stars on such highly eccentric orbits would get tidally disrupted by the SgrA* black hole, 
possibly producing both S-stars near the black hole and high-velocity stars in the Galactic halo. 
Subject headings: Stellar dynamics — galaxy: center — methods; N-body simulations 



1. INTRODUCTION 

Most well-observed galaxies with stellar bulges show 
evidenc e for supermassiv e black holes (SMBHs) in their 
centers dKormendvl |2004|) . An important feature of stellar 
dynamics in galactic nuclei is that for orbits contained well 
within the SMBH's radius of influence, the precession time 
iprcc is many orders of magnitude larger than their period P. ' 
Therefore, for a study of long-term (secular) gravitational 
dynamics, it is useful to think of the stars not as particles but 
as massive elliptical wires stretched along the stars' eccentric 
orbits (Rauch & Tremaine 1996). Gravitational torques 
between the wires can become the dominant mode of orbital 
eccentricity evolution, as is seen in the case of resonant re- 
laxation ([Ra uch & Tremaine 1996; Giir kan & Hopmanll2007l ; 
iTouma et al.ii2009k lEilon et al.ii2008i) . Secular effe cts have a 
significant impact on the stability of stel l ar cusps (|Tremainel 



2005 



2008) 



iPolyachenko & ShukhmanI I198U iPolyachenko et alJ 
on gravitational wave emission from in-spiraling 
compact objects (jHopman & Alexande ^2006alfbh . on tidal 
disruption of starsT Rauch & Ingallsu 19981) and, as we show 
here, on the origin of the S-stars in the Galactic center. 

In this Letter we describe a mechanism for generating high- 
eccentricity stars via a fast secular instability of an eccentric 
stellar disk (©. We use iV-body simulations to investigate 
the non-linear development of the instability (^and ® and 
to explore its dependence on the initial eccentricity and mass 
of the disk. We then show that the instability may be relevant 
for young stars in the Galactic center, and may provide an 
interesting channel for formation of the S-stars and of high- 
velocity stars in the Galactic halo (©. We critically discuss 
our results in ^ 

2. THE ECCENTRICITY INSTABILITY 

Consider a thin stellar disk of mass Mdisk around a SMBH 
of mass M, which is embedded in a power-law stellar cusp of 
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' This is in contrast to the rapidly precessing stellar orbits typical for glob- 
ular clusters or galaxies. 



mass Mcusp- Our two main assumptions are that (1) -Mdisk ^ 
^^cusp, so that precession of individual orbits is driven by the 
cusp, and (2) initially the stellar eccentricity vectors 
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(1) 



are aligned and similar in magnitude so that the disk pre- 
cesses coherently as a whole (for a typical power-law cusp, 
the orbital precession rate depends only weakly on the semi- 
major axis; see below). Here |e| corresponds to the magnitude 
of the eccentricity of the orbit and r and v are the position 
and velocity vectors respectively. Both assumptions are re- 
alistic for a disk of young stars formed in a galactic center 
via the gravitational capture of part of an infalling molecu- 
lar cloud; o bservations rev eal a massive spherical stellar cusp 
jGenzel et a l. 2003; Scho del et alj|2007l) . while hydrodynam- 
ical simulations tentatively show similar initial orbital eccen- 
tricities of young stars ^ ([S anders 1998; Alexander et al. 200^ 
iBonnell & Rici2008l;lHob bs & Nayakshin 2009). 

The orbits precess in the direction opposite to the orbital 
rotation of the stars. The precession time-scale scales as 
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where P{a) = 27r\/ a^/GM, is the period of a star with 
semi-major axis a, N{< a) is the number of stars in the cusp 
within a, m is the individual mass of the stars, and 7 is the 
power-law index for the spac e-density p(r ) oc r For a 
power-law cusp with 7 = 3/2 (lYoung|1980h . tprcc is constant 
for all a. 

Observations of the Galactic center indicate that there is a 
stellar cusp in the inner parsec, with published measurements 



^ IBonnell & Ricd 420081) make two large simulations, investigating star- 
formation in the Galactic center In their Run 1 they find 0.6 < e < 0.7, 
while in Run 2 the inner edge of the disk is circularized, and hence the stel- 
lar eccentricities range between in the inner disk and ~ 0.6 in the outer. 
While the precise nature of the circulaiization, or its absence, has not been 
investigated, it is clear that it depends on the initial conditions. 
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of the slope not significantly differe nt from 7 = 3/2, e.g. 
7 = 1.4 ± 0.1 (iGenzeletaLl lliol and 7 = 1.19 ± 0.05 
jSchodel et al.ll2007h . It is important to stress that these ob- 
servations only show a small and biased subset of the stellar 
content in that region. They are restricted to massive young 
stars, that have not yet relaxed, and to giants, which have 
relaxed, but are affected by hydr odynamical collisions with 
other stars (e.g. I Alexanden 1 1 9991) . These two stellar types 
are not expected to dominate the density. Instead, theoretical 
models accounting for mass-segregation show that the den- 
sity at the regions of interest is mostly determined by stel- 
lar black holes and main sequence s tars, with cusp values 
betwe en 7ms = 1-4 and 7bh = 2.0 ( [Alexander & HopmanI 
I2OO8I) . We conclude that both theory and observations yield 
values of 7 close to 3/2 and the dependence of tp^cc on a is 
weak. 

A key element of the eccentric instability is that /(e) in 
Equation dill is an increa sing function of eccentricity (e.g. 
iBinnev & Tremain^ l2008h . so that an orbit that is slightly 
more eccentric than average will lag behind in precession. 
Such a star will feel a strong, coherent torque from the other 
stars in the disk, in the opposite direction of its angular mo- 
mentum vector As a result, its angular momentum decreases 
in magnitude, causing its eccentricity to increase even fur- 
ther In this way, very high eccentricities can potentially be 
achieved. Conversely, if a stellar orbit is slightly less eccen- 
tric than the average, it will have a higher precession rate, thus 
experiencing a torque which acts to decrease its eccentricity 
further 

We thus find that the initial eccentricity vector distribution 
is unstable, due to a combination of retrograde precession and 
df /de > 0, and that all eccentricities are driven away from 
the initial eccentricity of the disk. It can be shown analyti- 
cally that the initial growth timescale of the eccentric instabil- 
ity scales as 
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(3) 

For the Galactic center this amounts to several precession 
time-scales. 

We now demonstrate the eccentricity instability using A^- 
body simulations and look at its impact on the young stars in 
the Galactic center. 

3. AT-BODY SIMULATIONS 

We summarize a newly-developed code which we have 
used in these simulations. Our in tegrator is based on a 4th- 
order Wisdom -Holman algorithm d Wisdom & Holmanll 19911: 
lYoshidalll990l) . where the forces on a particle are split into 
a dominant Keplerian force and perturb ation forces. We use 
Kepler's equation to integrate the orbit (lDanbvl[l992l) . keep- 
ing the position of the central object fixed. As we use adap- 
tive time-stepping to deal with close encounters, the algorithm 
loses its symplecticity and is not strictly energy conserving. 
We implement time-symmetr y to re duce the energy error and 
the K2 kernel from Deh nenI (1200 It) for gravitational soften- 
ing. The overall fractional energy error for these simulations 
is typically ^ lO^'* or less. We describe the code in detail in 
an upcoming paper 

3.1. Initial Conditions 

We have three main components in our simulations. (1) 
An eccentric disk with surface density profile S oc (cf. 



iPaumard et al.l2006l) . consisting of 100 equal-mass stars, with 
semi-major axes 0.05 pc < a < 0.5 pc. We vary the masses of 
the stars for different simulations. We initialize the eccentric- 
ity vectors of all stars pointing in roughly the same direction, 
though they are scattered in their orbital phase. Initial inclina- 
tions are dist ributed between h /r ~ ±0.01. (2) A SMBH of 
4 X 10^ Mq (iGhez et al.ll2008l) . and (3) a smooth stellar cusp 
with density profile i ndex 7 = 1.5 and a mass of 0.5 x l O^Mp, 
within 1 pc (Schod eretal]l2007HG"enzel et al.ll2003h . The 
stars within the disk precess due to the cusp and experience 
a prograde general relativistic precession which amounts to 
apsidal precession of the orbit by an angle 



6ttGM, 



a{l 



(4) 



per orbit. Precession due to self-gravity of the disk is minimal. 

To verify our results, we also run simulations with twice 
the number of half-as-massive stars and various softening pa- 
rameters. We d iscuss the use of different values of 7 (e.g. 
iBahcall & Wo"iiri976 ) in gg] 

4. RESULTS 

With a mass of 8 x IO'^Mq and an initial eccentricity of 
0.6, we let the disk of stars evolve for several precession times 
(ipioc ^ 0.6 - 0.7 Myrs or - 1200 orbits at a = 0.05 pc). In 
Figure ([T]) we show a selection of stars from this run which 
best demonstrates the dispersion of the disk's initial eccen- 
tricity, while Figure (|2|i follows the evolution of the eccentric- 
ity vectors (this time for A/disk = 4 x IO'^A/q). The circles 
in the latter plot indicate the stars experiencing the highest 
torque (r) and hence angular momentum changes. Note that 
the highest eccentricity orbits are precessing behind the bulk 
of the stars as expected. We find a monotonically decreasing 
relationship between the rate of change of angular momen- 
tum, i.e. torque on each orbit, and semi-major axis a. 




Fig. 1. — Eccentricity evolution of stars in an eccentric disk (Afjisk = 
8 x IO'^A/q) as a function of time. At is elapsed time expressed in units of 
the initial period of the innemiost orbit To (tproc ~ 1200 To). The smoothly 
varying behavior is characteristic of the secular nature of the instability. 

The effects of the instability on the stellar eccentricity 
distribution are evident after ^0.5 — 1 tprec and the strongest 
torques are experienced over a few precession times before 
the disk loses coherence. 
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Fig. 2. — Evolution of the eccentricity vectors of all the stars in the disk 
(Af(jisk = 4 X 10"^ Mq) for different times. The outermost circle indicates 
where |e| = 1; the inner circle indicates the initial value of the eccentricities 

e| = 0.6. An arrow points to the initial positions. The eccentricity vectors 
spread out as the stars complete more and more orbits, and precess with ret- 
rograde motion. Those experiencing the highest torques (t) are over-plotted 
with small circles. As the eccentric disk instabihty predicts, the slowest pre- 
cessing orbits have the highest eccentricities. In addition, orbits with low 
eccentricity clump together at late times. We believe non-linear effects are 
responsible for this feature. 



More puzzling in terms of origin is the popu lation of B-stars 
known as the 'S- stars' (Schodel et al. 2002; Ghez et alJl2005l: 
lEisenhauer et al.ll2005c iGillessen et al.ii2009i) . At distances 
of 0.003 — 0.03 pc from the SMBH, in-situ formation seems 
highly improbable due to the tidal field of SgrA* (iLevinI 
l2007h . In addition, their young age (-^ 20—100 Myr) imposes 
a tight constraint on formation scenarios in that they cannot 
have travelled far from their place of birth. One possibility 
is the form ation of the S-stars from the disru p tion of stel- 
lar bin aries ( lHillslll98 8; Gould & O uiiienll2003h . iPerets et alj 
(l2007h propose the surrounding stellar bulge as the source of 
the binaries, while iLockmann. Baumgardt. & Kroupal (l2008h 
(LBK) invoke two stellar disks in the Galactic center to ex- 
plain the origin. In their simulations the stars experienced 
Kozai-type torques which induced strong changes in their in- 
clination and eccentricity. The eccentricity of many of the 
stars became greater than 0.9, and LBK argued that if these 
stars were binaries they would be tidally disrupted by the 
SMBH, thus producing the S-stars and high-velocity stars via 
Hills' (119881) mechanism. 

However, the gravitational influence of a stellar cusp, 
which was neglected by LBK, strongly suppresses Kozai- 
type dynamics. This was recently discussed in lChangl(l2009l) . 
Chang's treatment was highly simplified; in particular, 
his expressions for Kozai evolution were valid only if the 
stellar orbits were inside the inner edges of both disks. We 
have performed LBK-type simulations and have confirmed 
Chang's general conclusion that a realistically massive stellar 
cusp entirely suppresses the production of high-eccentricity 
stars via a Kozai-type mechanism. 



We explore the dependence of the final range of eccentrici- 
ties on the mass and initial eccentricity of the disk. To this 
end, we run a set of simulations with a disk of initial eccen- 
tricity 0.6 for various disk masses (4, 6, 8, 10 x 10"^ A/0), and 
a set with iV/disk = 10"* M0 for various initial eccentricities 
(0. 1, 0.4, 0.5, 0.6, 0.7, 0.9). The results are presented in Table 
([T]). Higher initial eccentricities and masses of the disk are 
correlated with higher resulting eccentricities. We also see a 
correlation between the highest eccentricities and inclinations 
reached, with the largest values occurring at the smallest radii. 

In contrast to the secular instability described in 'Tremain^ 
(|2005), relativistic effects in this case do not promote stability 
once the relativistic prograde precession rate exceeds that of 
the retrograde precession due to the cusp. The highest eccen- 
tricity orbit precesses with prograde motion at a higher rate 
than the average and again experiences a torque in the oppo- 
site direction of its angular momentum. 

5. APPLICATION TO S-STARS 

T he region surroundin g the SMBH in the Galactic cen- 
ter (IGillessen et al.l 120091) is host to several interesting pop- 
ulations of young stars, which include one well-defined 
clockwise disk of O and Wolf-Rayet (WR) stars, and a 
group of counter-clockwise moving O and WR stars, pos- 
sibly a dissolving second disk dLev in & Beloborodov 20031 
Navakshin & Cuadra"2005VGenzel et al. 2003i:iPaumard et all 
2006; Lu et al. 2008; Bartko et al. 2008:^. These O and WR 
stars are located at projected distances of 0.05 pc to 0.5 pc, 
and were most likely formed as a result of gravitational f rag- 
mentation of a gaseous disk (iLevin & Beloborodovl2003l) . 



Table 1 

Simulation Parameters and Statistics of S-stars 





edisk 




ihigh'' 


Cbin 


rt 








am 


rp < rt'^ 


[Me] 






[rads] 


[AU] 


[pc] 








[pc] 


[%] 


4 X 10^ 


0. 4 


0.8755 


0.4662 


0.1 


3.57 


X 


10 


-5 


0.0026 













1 


3.57 


X 


10 


~4 


0.026 







0.6 


0.9691 


0.7960 


0.1 


3.57 


X 


10 


-5 


0.0026 













1 


3.57 


X 


10 


-4 


0.026 





6 X 10^ 


0.6 


0.9981 


1.824 


0.1 


3.57 


X 


10 


-5 


0.0026 


4 










1 


3.57 


X 


10 


-4 


0.026 
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8 X 10^ 


0. 4 


0.8796 


0.6643 


0.1 


3.57 


X 


10 


-5 


0.0026 













1 


3.57 


X 


10 


-4 


0.026 







0.6 


0.9993 


2.170 


0.1 


3.57 


X 


10 


~5 


0.0026 


3 
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3.57 


X 


10 


-4 


0.026 
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1 X 10** 


0.1 


0.5115 


0.2935 


0.1 


3.57 


X 


10 


~5 


0.0026 













1 


3.57 


X 


10 


-4 


0.026 







0.4 


0.9580 


1.440 


0.1 


3.57 


X 


10 


-5 


0.0026 


1 










1 


3.57 


X 


10 


-4 


0.026 
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0.5 


0.9773 


1.348 


0.1 


3.57 


X 


10 


~5 


0.0026 













1 


3.57 


X 


10 


-4 


0.026 


3 




0.6 


0.99991 


2.328 


0.1 


3.57 


X 


10 


-5 


0.0026 


14 










1 


3.57 


X 


10 


-4 


0.026 


14 




0.7 


0.99998 


3.0312 


0.1 


3.57 


X 


10 


-5 


0.0026 


32 










1 


3.57 


X 


10 


-4 


0.026 


32 




0.9 


1.0 


3.125 


0.1 


3.57 


X 


10 


-5 


0.0026 


58 










1 


3.57 


X 


10 


-4 


0.026 


59 



" Mean of highest five eccentricities reached during simulation 
^ Mean of highest five inclinations (in radians) 
Percentage of stars that enter the tidal radius 
Note that disks with A/disk 6 X 10^ Mq and initial eccentricity 
edisk ^ 0.6, produce very high eccentricities, possibly leading to tidal 
disruptions of binaries. 
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In contrast, the eccentricity instability is capable of driving 
several stars to near-radial orbits, as is evident in Figure 
([l]). For capture, the pericenter Tp of the binary's orbit must 

come within the tidal radius = (2A/,/r7ibin)^^'^abin- This 
results in a captured star with semi-major axis a, that scales 
as a, ^ (A/,/mbin)^^^abin, where Obin is the semi-major 
axis of the binary itself, TObin is its combined mass and the 
eccentricity of the captured star can be approximated as 

l-e~(mbi„/Af.)'/^- 

We summarize our results for varying masses and initial 
eccentricities of the disk in Table ([T]i- All simulations have 
evolved for ~ 7 Myr. For a given semi-major axis of a binary 
system, we calculate the tidal radius at which the binary will 
be disrupted and finally the percentage of stars that come 
within this radius. It is clear that for many parameters of 
our disk, stars have Vp < rt- For example, a disk with an 
initial eccentricity of 0.6 and mass of W^Mq, 14 % of stars 
in the disk will enter the tidal radius. Thus, if some of these 
stars are binaries. Hills' m echanism for producing the S- and 
high-velocity stars (iBrow n et al. 2006) could be at work. An 
important element in this scenari o is the bina ry survival rate 
at the radii of the disk; however, iPeretsI (l2009k) show that for 
a = 0.1 AU the evaporation rate is > lO'^ yr. 

The scalar resonant relaxation time in the Galactic center 
is everywhere larger than the age of the disks. It follows 
that since the post-capture eccentricity of the stars is of or- 
der 1 — e w 10"^, most of the S-stars cannot have formed 
in the current disk. The lifetime of many S-stars could be of 
the order of 100 Myr, and thus earUer 10-100 Myr) star- 
burst episodes could have contributed to their present popula- 
tion. In that case, the scalar re laxation time is short enough 
to randomize the e ccentricities dHopman & Alexanded2006al : 
iPerets et"ani2009h . We note that the scenario presented by 
LBK has the same limitation. 



5.1. Bimodal Eccentricity Distribution in Disk 
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Fig. 3. — (a) Observed eccentricity distribution of clock-wise stellar disk 
in the Galactic center, reproduced with permission from Bartko et al. (2008). 
(b) Simulated distribution for a disk with an initial eccentricity of 0.6 and 
A^disk = 10^ after 5.7 Myr. Both the observed and simulated eccentricity 
distributions exhibit bimodality and are in good agreement. 

In Figure (O we compare the observ ed stellar eccentric- 
ity distribution of the clock-wise disk dBartko et alj l2008l) 
with the simulated post-instability eccentricity distribution 
for all stars with an initial disk eccentricity of 0.6 (c.f. 
iBonnell & Ricell2008l) . We draw attention to the qualitative 
similarities in the two plots, in particular the double-peaked 



profile. This bimodality is a natural consequence of the eccen- 
tric disk instability as all eccentricities are driven away from 
the mean. Without the instability it is unclear how to generate 
such a distribution. 

6. CRITICAL DISCUSSION 

So far, we have described the eccentricity instability for 
an idealized initial disk configuration. We now demonstrate 
the effectiveness of the instability for a broader scope of 
parameters and extend our disk model to include (1) a 
significant range in both magnitude and direction of the 
stellar eccentricity vectors, (2) heavily inclined stellar orbits, 
and (3) different values of 7, the power-law index of the cusp, 
which determines the precession rate within the disk. We 
discuss each of these in turn. 

(1) We simulate disks with initial eccentricity vectors 
scattered between a range of opening angles 9 E [0, 27r]. 
As anticipated, significant spreading of the vector decreases 
the effectiveness of the instability. However, we find suf- 
ficiently high eccentricity orbits to generate S-stars up to 
9 ~ TT radians, which suggests that the instability remains 
important for substantially less favorable circumstances. 
Next we vary the magnitudes of the eccentricities throughout 
the disk, ba sing our initial conditions on Run 1 and 2 of 
IBonnell & R ice (2008). Run 1 (0.6 < e < 0.7) is in excellent 
agreement with our previous results. The instability in Run 
2 (0 < e < 0.6) however, is not efficient at generating the 
highest e orbits and hence the S-stars. These conclusions are 
transparent; the stars which experience the greatest torque are 
those at small a (see © and if initially circular, need a greater 
torque to produce the highest e orbits. We find that the most 
important condition, for Mdisk = 10^, is for innermost orbits 
to have e > 0.4. Future work on gas dynamics of eccentric 
disks will tell whether this presents a serious problem for our 
scenario. 

(2) We study disks with stellar inclinations h/r in the range 
[—0.3, 0.3] and find that the instability is robust in all cases. 
However, above the limit [—0.2,0.2], the instability appears 
ineffective in generating the highest e (> 1 — 10^'^) orbits. 

(3) We vary the power-law slope of the cusp. For 7 < 3/2, 
due to the shallowness of the cusp, the innermost orbits 
lag behind in precession and are more susceptible to being 
pushed to high e orbits. While this improves the ability of the 
instability to produce S-stars, the instability is not as effective 
at generating low e orbits. Conversely, an increase in 7 
(> 3/2) results in the innermost orbits precessing faster than 
the bulk of the stars and hence suppresses the generation of 
the highest e orbits and lengthens the time-scale on which it 
occurs. We conclude that the production of S-stars is unlikely 
for a cusp with 7 > 7/4. In all cases examined (1 < 7 < 2), 
the instability proves robust and, as emphasized in ^ there 
is both observational and theoretical evidence that the slope 
of the cusp in the Galactic center is sufficiently close to 3/2 
for the instability to be highly effective. 

The eccentricity instability described in this work is rele- 
vant for an eccentric disk embed ded in a stellar cu sp. We note 
that the eccentric disk in M31 (lTremainelll995h is not sub- 
ject to our instability, since it is not embedded in any visible 
cusp. Indeed we find that such a disk remains coherent over 
several precession times, although the innermost stars can un- 
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dergo significant an gular momentum changes as described by 
ICuadra et alJ (l2008h . The situation is different in the Galactic 
center, where massive stars were likely born from an eccen- 
tric gaseous disk, and where there is strong obse rvational ev- 
idence for a heavy spherical stellar cusp (IGenzel et alJl2003l : 
ISchodel eFal. 2007) . We have shown that the instability may 
have strongly influenced the development of the eccentricities 
of the massive stars observed, and in particular, can account 
for the bimodality in their e ccentr icity distribution. Further- 
more, combined with Hills' (1 19881) mechanism, the eccentric- 
ity instability may have played a role in the formation of the 



S-stars and of the high-velocity stars in the Galactic halo. 
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